# Subset only
# This file re-runs HPBM results
# and our replication of HPBM at agency and county level
# and stores all of the results
library(plm) #for panel data
library(lfe) #for felm function
library(rlang)
library(tidyverse)
library(broom)
library(foreign)
library(purrr)
library(texreg)
library(foreign)
source("lib/model-functions-HPBM.R") #these have all the different models we need to run
options(scipen=0, digits=7) #to standardize number of digits when exporting data. all data here is run using these options.

year_subset <- TRUE
city_outliers_eliminated <- TRUE

# load data
load('data/city-aggregate.RData')
load('data/county-aggregate.RData')
load(file='data/harris.RData') #Harris data for exact replication
harris <- df
rm(df)

## load controlled data
#load('data/city-aggregate-controlled.RData')
#load('data/county-aggregate-controlled.RData')

# restrict data
if(year_subset){
  agency <- agency[agency$year >=2010 & agency$year<=2012,]
  county <- county[county$year >=2010 & county$year <=2012,]
  harris <- harris[harris$year>=2010 & harris$year<=2012,]
}

if(city_outliers_eliminated){
  agency <- agency %>% filter(allcrime_rate < 1000000)
}

# create panel data
agency <- pdata.frame(agency,index=c("g","year"))
county <- pdata.frame(county,index=c("g","year"))
harris <- pdata.frame(harris, index = c("fips", "year")) 

#------------------------------------------------------------#
## Harris tables for supplementary materials ----
#------------------------------------------------------------#

## Harris Replication using our Agency level data -----
homicide.items.agency.Harris <- run.IV.agency.Harris(Murder_rate, agency, 'items')
homicide.value.agency.Harris <- run.IV.agency.Harris(Murder_rate, agency, 'value')
robbery.items.agency.Harris <- run.IV.agency.Harris(Robbery_rate, agency, 'items')
robbery.value.agency.Harris <- run.IV.agency.Harris(Robbery_rate, agency, 'value')
assault.items.agency.Harris <- run.IV.agency.Harris(Assault_rate, agency, 'items')
assault.value.agency.Harris <- run.IV.agency.Harris(Assault_rate, agency, 'value')
vehicle.items.agency.Harris <- run.IV.agency.Harris(MV_theft_rate, agency, 'items')
vehicle.value.agency.Harris <- run.IV.agency.Harris(MV_theft_rate, agency, 'value')


homicide.items.agency.Harris <- run.IV.agency.Harris(Murder_rate, agency, 'items')
homicide.value.agency.Harris <- run.IV.agency.Harris(Murder_rate, agency, 'value')
robbery.items.agency.Harris <- run.IV.agency.Harris(Robbery_rate, agency, 'items')
robbery.value.agency.Harris <- run.IV.agency.Harris(Robbery_rate, agency, 'value')

# reduced form at agency level
homicide.red.agency.items.harris <- reduced.agency.Harris(Murder_rate, agency, 'items')
homicide.red.agency.value.harris <- reduced.agency.Harris(Murder_rate, agency, 'value')
robbery.red.agency.items.harris <- reduced.agency.Harris(Robbery_rate, agency, 'items')
robbery.red.agency.value.harris <- reduced.agency.Harris(Robbery_rate, agency, 'value')
assault.red.agency.items.harris <- reduced.agency.Harris(Assault_rate, agency, 'items')
assault.red.agency.value.harris <- reduced.agency.Harris(Assault_rate, agency, 'value')
vehicle.red.agency.items.harris <- reduced.agency.Harris(MV_theft_rate, agency, 'items')
vehicle.red.agency.value.harris <- reduced.agency.Harris(MV_theft_rate, agency, 'value')

## Harris County Replication -------
homicide.county.items.harris <- run.IV.county.Harris(MURDER_countysumrate, county, 'items')
homicide.county.value.harris <- run.IV.county.Harris(MURDER_countysumrate, county, 'value')
robbery.county.items.harris <- run.IV.county.Harris(ROBBERY_countysumrate, county, 'items')
robbery.county.value.harris <- run.IV.county.Harris(ROBBERY_countysumrate, county, 'value')
assault.county.items.harris <- run.IV.county.Harris(AGASSLT_countysumrate, county, 'items')
assault.county.value.harris <- run.IV.county.Harris(AGASSLT_countysumrate, county, 'value')
vehicle.county.items.harris <- run.IV.county.Harris(MVTHEFT_countysumrate, county, 'items')
vehicle.county.value.harris <- run.IV.county.Harris(MVTHEFT_countysumrate, county, 'value')


# reduced form at county level
homicide.red.county.items.harris <- reduced.county.Harris(MURDER_countysumrate, county, 'items')
homicide.red.county.value.harris <- reduced.county.Harris(MURDER_countysumrate, county, 'value')
robbery.red.county.items.harris <- reduced.county.Harris(ROBBERY_countysumrate, county, 'items')
robbery.red.county.value.harris <- reduced.county.Harris(ROBBERY_countysumrate, county, 'value')
assault.red.county.items.harris <- reduced.county.Harris(AGASSLT_countysumrate, county, 'items')
assault.red.county.value.harris <- reduced.county.Harris(AGASSLT_countysumrate, county, 'value')
vehicle.red.county.items.harris <- reduced.county.Harris(MVTHEFT_countysumrate, county, 'items')
vehicle.red.county.value.harris <- reduced.county.Harris(MVTHEFT_countysumrate, county, 'value')


## Exact replication of Harris Results -------
homicide.items.harris <- run.IV.Harris(Amurders, harris, 'items')
homicide.value.harris <- run.IV.Harris(Amurders, harris, 'value')
robbery.items.harris <- run.IV.Harris(Arobbery, harris, 'items')
robbery.value.harris <- run.IV.Harris(Arobbery, harris, 'value')
assault.items.harris <- run.IV.Harris(Aassault, harris, 'items')
assault.value.harris <- run.IV.Harris(Aassault, harris, 'value')
vehicle.items.harris <- run.IV.Harris(Avehicle, harris, 'items')
vehicle.value.harris <- run.IV.Harris(Avehicle, harris, 'value')

homicide.red.items.harris <- reduced.Harris(Amurders, harris, 'items')
homicide.red.value.harris <- reduced.Harris(Amurders, harris, 'value')
robbery.red.items.harris <- reduced.Harris(Arobbery, harris, 'items')
robbery.red.value.harris <- reduced.Harris(Arobbery, harris, 'value')
assault.red.items.harris <- reduced.Harris(Aassault, harris, 'items')
assault.red.value.harris <- reduced.Harris(Aassault, harris, 'value')
vehicle.red.items.harris <- reduced.Harris(Avehicle, harris, 'items')
vehicle.red.value.harris <- reduced.Harris(Avehicle, harris, 'value')

#---------------------------------------------------------#
## First stage --------
#---------------------------------------------------------#
harris_agency_items_first <- summary(homicide.items.agency.Harris$stage1)$coefficients
harris_agency_value_first <- summary(homicide.value.agency.Harris$stage1)$coefficients
harris_county_items_first <- summary(homicide.county.items.harris$stage1)$coefficients
harris_county_value_first <- summary(homicide.county.value.harris$stage1)$coefficients
harris_items_first <- summary(homicide.items.harris$stage1)$coefficients
harris_value_first <- summary(homicide.value.harris$stage1)$coefficients
first_list <- list(harris_agency_items_first,
                   harris_agency_value_first,
                   harris_county_items_first,
                   harris_county_value_first,
                   harris_items_first,
                   harris_value_first)
HPBM_first_stage_results <- tibble(first_stage = map(first_list, tidy))

# Save all results ------
agency_items_list <- list(homicide.items.agency.Harris,
                          robbery.items.agency.Harris,
                          assault.items.agency.Harris,
                          vehicle.items.agency.Harris)
agency_value_list <- list(homicide.value.agency.Harris,
                          robbery.value.agency.Harris,
                          assault.value.agency.Harris,
                          vehicle.value.agency.Harris)

county_items_list <- list(homicide.county.items.harris,
                          robbery.county.items.harris,
                          assault.county.items.harris,
                          vehicle.county.items.harris)
county_value_list <- list(homicide.county.value.harris,
                          robbery.county.value.harris,
                          assault.county.value.harris,
                          vehicle.county.value.harris)


HPBM_agency_items_results <- tibble(map(agency_items_list, tidy))
HPBM_agency_value_results <- tibble(map(agency_value_list, tidy))
HPBM_county_items_results <- tibble(map(county_items_list, tidy))
HPBM_county_value_results <- tibble(map(county_value_list, tidy))

agency_red_items_list <- list(homicide.red.agency.items.harris,
                              robbery.red.agency.items.harris,
                              assault.red.agency.items.harris,
                              vehicle.red.agency.items.harris)
agency_red_value_list <- list(homicide.red.agency.value.harris,
                              robbery.red.agency.value.harris,
                              assault.red.agency.value.harris,
                              vehicle.red.agency.value.harris)

county_red_items_list <- list(homicide.red.county.items.harris,
                              robbery.red.county.items.harris,
                              assault.red.county.items.harris,
                              vehicle.red.county.items.harris)
county_red_value_list <- list(homicide.red.county.value.harris,
                              robbery.red.county.value.harris,
                              assault.red.county.value.harris,
                              vehicle.red.county.value.harris)


HPBM_agency_red_items_results <- tibble(map(agency_red_items_list, tidy))
HPBM_agency_red_value_results <- tibble(map(agency_red_value_list, tidy))
HPBM_county_red_items_results <- tibble(map(county_red_items_list, tidy))
HPBM_county_red_value_results <- tibble(map(county_red_value_list, tidy))


harris_items_list <- list(homicide.items.harris,
                          robbery.items.harris,
                          assault.items.harris,
                          vehicle.items.harris)
harris_value_list <- list(homicide.value.harris,
                          robbery.value.harris,
                          assault.value.harris,
                          vehicle.value.harris)

harris_red_items_list <- list(homicide.red.items.harris,
                              robbery.red.items.harris,
                              assault.red.items.harris,
                              vehicle.red.items.harris)
harris_red_value_list <- list(homicide.red.value.harris,
                              robbery.red.value.harris,
                              assault.red.value.harris,
                              vehicle.red.value.harris)

HPBM_items_results <- tibble(map(harris_items_list,tidy))
HPBM_value_results <- tibble(map(harris_value_list,tidy))
HPBM_red_items_results <- tibble(map(harris_red_items_list,tidy))
HPBM_red_value_results <- tibble(map(harris_red_value_list,tidy))

# save
save(HPBM_agency_items_results,
     HPBM_agency_value_results,
     HPBM_county_items_results,
     HPBM_county_value_results,
     HPBM_agency_red_items_results,
     HPBM_agency_red_value_results,
     HPBM_county_red_items_results,
     HPBM_county_red_value_results,
     HPBM_items_results,
     HPBM_value_results,
     HPBM_red_items_results,
     HPBM_red_value_results,
     HPBM_first_stage_results,
     file = 'data/results-HPBM.RData')

# Save all results for texreg
HPBM_agency_items_results <- map(agency_items_list, texreg::extract)
HPBM_agency_value_results <- map(agency_value_list, texreg::extract)
HPBM_county_items_results <- map(county_items_list, texreg::extract)
HPBM_county_value_results <- map(county_value_list, texreg::extract)

HPBM_agency_red_items_results <- map(agency_red_items_list, texreg::extract)
HPBM_agency_red_value_results <- map(agency_red_value_list, texreg::extract)
HPBM_county_red_items_results <- map(county_red_items_list, texreg::extract)
HPBM_county_red_value_results <- map(county_red_value_list, texreg::extract)

HPBM_items_results <- map(harris_items_list, texreg::extract)
HPBM_value_results <- map(harris_value_list, texreg::extract)
HPBM_red_items_results <- map(harris_red_items_list, texreg::extract)
HPBM_red_value_results <- map(harris_red_value_list, texreg::extract)

harris_agency_items_first <- homicide.items.agency.Harris$stage1
harris_agency_value_first <- homicide.value.agency.Harris$stage1
harris_county_items_first <- homicide.county.items.harris$stage1
harris_county_value_first <- homicide.county.value.harris$stage1
harris_items_first <- homicide.items.harris$stage1
harris_value_first <- homicide.value.harris$stage1
first_list <- list(harris_agency_items_first,
                   harris_agency_value_first,
                   harris_county_items_first,
                   harris_county_value_first,
                   harris_items_first,
                   harris_value_first)
HPBM_first_stage_results <- map(first_list, texreg::extract)


# save subsetted results
save(HPBM_agency_items_results,
     HPBM_agency_value_results,
     HPBM_county_items_results,
     HPBM_county_value_results,
     HPBM_agency_red_items_results,
     HPBM_agency_red_value_results,
     HPBM_county_red_items_results,
     HPBM_county_red_value_results,
     HPBM_items_results,
     HPBM_value_results,
     HPBM_red_items_results,
     HPBM_red_value_results,
     HPBM_first_stage_results,
     file = 'data/results-HPBM-for-tables-subset.RData')

# save results with outliers included
save(HPBM_agency_items_results,
     HPBM_agency_value_results,
     HPBM_county_items_results,
     HPBM_county_value_results,
     HPBM_agency_red_items_results,
     HPBM_agency_red_value_results,
     HPBM_county_red_items_results,
     HPBM_county_red_value_results,
     HPBM_items_results,
     HPBM_value_results,
     HPBM_red_items_results,
     HPBM_red_value_results,
     HPBM_first_stage_results,
     file = 'data/results-HPBM-for-tables-outliers.RData')

# save results for controlled items
save(HPBM_agency_items_results,
     HPBM_agency_value_results,
     HPBM_county_items_results,
     HPBM_county_value_results,
     HPBM_agency_red_items_results,
     HPBM_agency_red_value_results,
     HPBM_county_red_items_results,
     HPBM_county_red_value_results,
     HPBM_items_results,
     HPBM_value_results,
     HPBM_red_items_results,
     HPBM_red_value_results,
     HPBM_first_stage_results,
     file = 'data/results-HPBM-for-tables-controlled.RData')



